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APPLICATION  OF  NEW  ALGORITHMS  FOR 
THE  SIMULATION  OF  MULTIDISCIPLINARY  PROBLEMS: 
FLUID,  STRUCTURE,  THERMAL  AND  CONTROL  COUPLING 


Rainald  Lohner  and  Chi  Yang 
School  of  Computational  Sciences 
George  Mason  University 
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Summary 

The  overall  objective  of  the  research  carried  out  over  the  last  three  yezirs  was  the 
development  of  new  algorithms  for  the  efficient  simulation  of  multidisciplinary  problems 
that  require  the  simulation  of  viscous,  conducting,  compressible  flows  around  or  within 
moving,  structurally  and  thermally  responding  structures.  The  development  was  based 
on  current  3-D  Euler/Navier-Stokes  capabilities,  and  encompassed  flow  solvers,  grid 
generation,  CFD-CSD-CTD  integration,  optimal  shape  design,  and  the  efficient  use  of 
emerging  super  computer  hardware.  The  research  carried  out  over  the  last  three  years 
significantly  advanced  the  state  of  the  art  in  this  area  of  CFD.  The  particular  topics 
are  treated  below  in  detail. 
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1.  FLOW  SOLVERS 

For  the  flow  solvers,  seven  major  developments  took  place  over  the  course  of  this  re¬ 
search  effort: 

-  Incompressible  solver  modules: 

-  New  implicit-advection  projection  solver; 

-  Pressure  boundary  conditions  for  internal  flows; 

-  Arbitrary  number  of  advective/diffusive  species; 

-  Inclusion  of  vorticity  confinement  to  track  vortices  over  large  distances; 

-  Combination  of  Baldwin-Lomax  and  Smagorinsky  turbulence  models; 

-  Compressible  solver  modules: 

-  Two  matrix-free  implicit  compressible  solvers; 

-  Arbitrary-gas  equation  of  state  lookup; 

-  Combustion  modeling. 

1.1  Incompressible  Flow  Solvers 

New  Implicit-Advection  Projection  Solver:  To  date,  the  advective  terms  present  in  the 
Navier-Stokes  equations  have  been  integrated  explicitly.  This  yields  highly  accurate 
results,  but  also  imposes  severe  restrictions  on  the  allowable  timestep.  For  grids  that 
exhibit  a  large  variation  in  element  size  the  timestep  imposed  by  the  smallest  elements 
may  be  orders  of  magnitude  smaller  than  the  timestep  required  to  obtain  time-accurate 
results.  This  implies  that  for  many  classes  of  problems,  e.g.  heat  release  with  buoyancy 
effects,  tens  of  thousands  of  timesteps  are  required  per  simulation,  rendering  the  scheme 
impractical.  This  prompted  the  quest  for  an  implicit  integration  of  the  advective  terms 
that  could  be  incorporated  easily  into  the  existing,  mature,  explicit-advection  projec¬ 
tion  scheme. 

A  clciss  of  implicit  projection  schemes  was  developed  that  allows  to  circumvent  the 
stringent  explicit  advection  timestep  limitation  without  incurring  a  major  code  change. 
The  core  solver  (explicit  advection  with  local  timestepping)  is  maintained  intact,  with 
the  exception  of  simple  source-terms  that  are  pointwise  and  therefore  themselves  well  for 
implicit  integration.  The  implicit  projection  schemes  were  implemented  into  FEFLO, 
a  general  pmpose  CFD  code  based  on  unstructured-grids  and  simple  finite  elements. 
The  implicit  projection  schemes  were  then  used  to  solve  transient  flows,  some  of  which 
had  thermal  buoyancy  effects.  The  results  show  considerable  savings  for  the  implicit 
schemes  as  compared  to  their  explicit  advection  counterparts. 

For  an  application  which  is  a  time-accurate  run,  the  new  implicit  scheme  runs  approx¬ 
imately  10  times  faster  than  the  explicit  scheme.  The  option  of  advecting  an  arbitrary 
nmnber  of  species  was  added  to  simulate  the  effects  of  contaminants,  combustion  and 
for  diagnostics  purposes. 

The  new  scheme  has  been  tested  for  a  number  of  benchmark  cases.  Some  of  the  test 
cases  are  listed  as  follows: 

-  VonKarman  vortex  street,  cooling  tower; 
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-  Natural  convection  of  air  in  a  square  cavity 

-  Cooling  tower; 

-  Steady  ship  wave  around  Wigley  hull. 

Pressure  Boundary  Conditions  For  Internal  flows:  For  external  aerodynamics,  the  im¬ 
position  of  pressure  at  outflow  boundaries  presents  no  further  difficulties.  The  situation 
is  very  different  for  internal  flows,  such  as  pipe  networks  or  arteries  in  the  human  body. 
Two  new  pressure  boundary  conditions  for  internal  flows  were  derived  and  implemented 
in  FEFLO. 

It  was  shown  in  the  Circle  of  Willis  in  the  brain  the  new  pressure  boundary  conditions 
were  the  only  way  to  simulate  with  some  degree  of  realism  the  complex  hemodynamic 
patterns  that  arise. 

Vorticitv  Confinement:  A  new,  general  vorticity  confinement  term  has  been  derived 
using  dimensional  analysis.  The  resulting  vorticity  confinement  term  is  a  function  of 
the  local  vorticity-based  Reynolds  number,  the  local  element  size,  the  vorticity  and  the 
gradient  of  the  vorticity.  The  vorticity  confinement  term  disappears  for  vanishing  mesh 
size,  and  is  applicable  to  unstructmed  grids  with  large  element  size  disparity. 

The  new  vorticity  confinement  term  was  implemented  into  FEFLO.  Different  vorticity 
confinement  options  have  been  tested  and  compared,  and  the  best  was  selected.  The 
new  term  has  been  found  to  be  successful  for  a  number  of  test  cases,  allowing  better 
definition  of  vortices  without  any  deleterious  effects  on  the  flow  field. 

The  verification  and  validation  are  carried  out  for  the  following  cases: 

-  von-Karman  vortex  street; 

-  Wing  tip  vortex  separation; 

-  Wing  vortex  interaction;  and 

-  Massively  separated  flow  over  a  wing. 

The  vortex  core  visualization  for  the  NACA0012  wing  showed  that  the  vortex  was  well 
captmed  till  eight  chord  lengths  downstream  of  the  trailing  edge  when  the  vorticity 
confinement  term  was  switched  on.  The  vortex  was  dissipated  after  only  one  chord 
length  without  vorticity  confinement.  The  vorticity  confinement  also  had  a  marked 
effect  on  the  strength  of  the  detached  vortex  for  the  delta  wing  case,  and  the  strength 
of  the  von-Karman  vortex  for  the  2-D  cylinder.  The  large  scale  demonstration  was 
performed  for  a  submarine  with  a  42-foot  diameter.  The  speed  was  7  knots  and  the 
Reynolds  munber  was  1.1x10®. 

Baldwin-Lomax-Smagorinskv  Model:  The  path  to  high-fidelity,  high  Reynolds-number 
simulations  always  leads  to  impossibly  high  numbers  of  gridpoints  if  all  temporal  and 
spatial  wavelengths  present  in  the  flow  are  resolved.  In  order  to  obtain  tractable  num¬ 
bers  of  gridpoints,  either  temporal  averaging  (via  Reynolds- Averaged  Navier-Stokes 
(RANS))  or  some  form  of  spatial  averaging  (via  Large  Eddy  Simulations  (LES))  is  re¬ 
quired.  A  relative  new  and  promising  approach  is  Detached  Eddy  Simluation  (DES), 
which  merges  the  best  of  RANS  and  LES.  RANS  is  used  near  the  wall  where  there  is 
no  separation.  LES  is  used  in  the  core  flow  where  it  works  well.  In  this  spirit,  the  two 
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simplest  algebraic  models:  Baldwin  Lomax  and  Smagorinsky  were  combined  in  such  a 
way  that  each  one  is  used  in  the  region  of  the  flow  it  was  devised  for.  The  Baldwin 
Lomax  (BL)  turbulence  model  is  used  close  to  walls,  where  it  is  known  to  give  very 
good  results.  Smagorinsky  (SMA)  is  used  for  regions  of  separation  and  core  flow.  In 
a  first  pass,  the  SMA  viscosity  is  computed  from  the  vorticity  or  second  invariant  of 
the  deformation  rate  tensor,  as  well  as  the  local  mesh  size.  In  a  second  pass,  the  BL 
viscosity  is  computed  as  usual,  i.e.  by  dividing  the  boundary  layer  into  an  inner  and 
outer  layer.  For  the  inner  layer,  the  BL  viscosity  is  kept.  In  the  outer  layer,  the  BL 
viscosity  is  used  until  it  falls  below  the  value  of  the  SMA  viscosity.  In  the  remainder 
of  the  flow  domain,  the  SMA  viscosity  is  used. 

1.2  Compressible  Flow  Solvers 

Two  Matrix-Free  Implicit  Compressible  Solvers:  If  we  consider  a  steady  flow  problem, 
the  fastest  solvers  to  date  include  multigrid  and  GMRES-LU-SGS  solvers.  For  super¬ 
sonic  flows,  space-marching  and  blocking  represent  an  interesting  alternative.  During 
the  past  three  years,  the  GMRES-LU-SGS  solvers  developed  in  prototype  versions  of 
FEFLO  were  migrated  into  the  production  code. 

The  classic  relaxation  schemes  to  solve  the  resulting  system  have  been  optimized  over 
the  years,  resulting  in  very  efficient  solvers.  The  combined  effect  of  these  simplifications 
is  a  family  of  schemes  that  are  matrix  free,  require  no  extra  storage  as  compared  to 
explicit  schemes,  and  per  relaxation  sweep  are  faster  than  conventional  expUcit  schemes. 

Abitrarv-Gas  Equation  of  State  Lookup:  In  order  to  be  able  to  model  effectively  mix¬ 
tures  of  gases  such  as  those  encountered  in  nozzle  exhaust  plumes,  a  conversion  code 
was  written  for  the  NASA-Glenn  chemical  equilibrium  program  CEA  of  Gordon  and 
McBride.  The  resulting  table  is  output  in  the  Los  Alamos  SESAME  format,  and  read 
in  at  the  beggining  of  the  run.  The  table  look-up  reqxiires  less  than  5%  of  the  total 
CPU  time  for  expUcit  solvers  on  RISC  machines,  but  has  not  yet  been  ported  to  vector 
machines. 

Gnmbnstion  Modeling:  During  past  three  years,  we  implemented  both  equilibrium  and 
finite  rate  models.  The  core  modules  follow  closely  the  techniques  of  CONCHAS  and 
KIVA.  The  stiff  ODE  integrator  decomposes  the  reactions  into  fast  and  slow  ones,  and 
applies  different  integration  techniques  depending  on  the  speed. 

2.  GRID  GENERATION 

In  the  area  of  grid  generation,  there  were  five  major  developments  that  took  place 
during  the  course  of  this  research  effort: 

-  Advances  in  the  gridding  of  surfaces  given  as  discrete  data; 

-  Improvements  in  RANS  gridding; 

-  Improvements  in  parallel  grid  generation. 

2.1  Discrete  Surface  Meshing 
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The  rapid,  user-friendly  definition  of  the  surfaces  defining  the  computational  domain 
has  been  an  important  goal  during  the  last  decade.  Surfaces  can  be  defined  either  an¬ 
alytically  (using  B-Splines,  NURBS,  Coon’s  patches,  etc.)  or  via  triangulations.  The 
latter  option  is  particularly  interesting  for  data  sets  stemming  from  remote  sensing  data 
(e.g.  geographical  data)  or  medical  imaging.  An  interesting  observation  made  over  the 
last  years  is  that  an  increasing  number  of  data  sets  used  to  define  the  geometry  of  CFD 
domains  is  given  in  the  form  of  triangulations,  even  though  the  CAD  data  is  available. 
The  reason  for  this  shift  in  data  type  is  that  a  watertight  triangulation  defines  in  a 
unique  way  the  domain  considered,  and  does  not  require  any  futher  geometric  cleanup 
operations.  This  is  not  the  case  with  native  CAD  datasets,  in  which  we  frequently 
encounter  very  large  numbers  of  patches,  overlapping  patches,  gaps,  and  other  geo¬ 
metric  pathologies  that  require  user  intervention.  These  developments  have  renewed 
the  interest  in  robust  surface  meshing  of  so-called  discrete  surfaces  (DS).  Of  the  many 
innovations  introduced  during  the  last  years,  we  mention; 

-  Automatic  preprocessing/improvement  of  the  DS; 

-  Introduction  of  a  visibility  horizon  filter  for  close  points/  sides; 

-  Strict  enforcement  of  continuous  topology; 

-  Improved  2D  cross-check;  and 

-  Adaptive  background  grid  element  size  definition. 

In  the  sequel,  we  expand  on  a  few  of  these. 

Visibility  Horizon:  The  advancing  front  method  adds  a  new  surface  triangle  by  remov¬ 
ing  a  side  from  the  active  front.  Among  the  decisions  required  is  whether  to  take  an 
existing  point  to  form  the  new  triangle,  or  to  introduce  a  new  point.  The  list  of  close 
(i.e.  possible)  points  is  obtained  from  a  proximity  search.  This  list  of  possible  close 
points  is  reduced  by  several  tests  (visibility,  angles,  etc.).  Perhaps  the  most  important 
validation  test  is  based  on  the  neighbom  to  neighbour  search  on  the  given  DS.  The 
starting  face  for  the  search  is  given  by  the  underlying  DS  face  at  the  midpoint  of  the 
side  being  removed  from  the  active  front.  The  direction  is  given  by  the  close  point.  Any 
close  point  that  can  not  be  reached  on  the  given  DS  using  the  neighbour  to  neighbour 
search  is  removed  from  the  fist.  A  similar  procedure  is  used  to  filter  close  sides,  which 
are  required  to  test  if  the  new  triangle  crosses  the  existing  active  front  of  sides. 

Continuous  Topology:  A  typical  neighbom:  to  neighbour  search  will  not  stop  at  internal 
DS  geometry  lines  (given  by  sharp  edges).  Therefore,  a  so-called  visibility  horizon  was 
introduced  for  the  neighbom  to  neighbom  search.  All  neighbom  to  neighbom  edges 
given  by  internal  geometry  lines  or  angles  beyond  a  certain  tolerance  are  marked.  In 
this  way,  the  neighbom  to  neighbour  search  can  recognize  them.  The  neighbom  to 
neighbom  search  stops  at  these  internal  geometry  fines.  The  close  point  is  marked  as 
mireachable  and  removed  from  the  fist  of  candidates. 

Adaptive  Background  Grid  Based  on  DS:  For  complex  geometries,  the  specification  of 
desired  element  size  can  be  a  tedious,  time-consuming  process.  Adaptive  background 
grids,  offer  the  possibility  to  reduce  drastically  the  required  level  of  human  input.  DS 
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offer,  by  their  way  of  defining  the  surface,  a  natural  way  to  refine  the  background  grid 
and  to  define  the  mesh  size  required  for  a  proper  definition  of  the  geometry. 

2.2  Improvements  in  RANS  Gridding 

The  generation  of  high-quality  grids  suitable  for  RANS  calculations  of  flows  in  and 
around  complex  geometries  continues  to  be  an  active  area  of  research.  The  two  most 
common  ways  of  generating  highly  stretched  grids  suitable  for  RANS  calculations  are 
the  advancing  layers  technique  (ALT)  and  the  directional  enrichment  technique  (DET). 
The  ALT  follows  the  spirit  of  the  advancing  front  technique  (AFT):  starting  from  the 
‘wetted  surface’,  add  thin  layers  of  elements  until  an  isotropic  mesh  is  achieved.  Prom 
this  point  onwards,  the  mesh  is  completed  with  the  AFT.  The  DET  is  an  extension  of 
the  Generalized  Delaunay  Triangulation  (GDT).  At  first,  an  isotropic  mesh  is  generated 
(either  via  AFT  or  GDT).  Then,  the  points  in  the  regions  where  stretched  elements 
are  to  be  generated  are  removed  (typically  by  collapsing  edges).  Thereafter,  points 
are  introduced  in  the  near-wall  region  so  as  to  obtain  the  desired  RANS  mesh.  The 
newly  introduced  points  are  reconnected  using  the  GDT.  This  latter  technique,  while 
general,  can  still  generate  bad  grids  for  complex  geometries.  A  recent  improvement  was 
obtained  by  considering  the  layer-number  of  the  points  when  performing  diagonal/face 
swapping  in  3-D.  The  quality  of  elements  that  cross  several  point-layers  is  reduced, 
leading  to  swaps  that  favour  elements  with  minimal  layer  jump. 

2.3  Parallel  Grid  Generation 

Over  the  last  decade,  major  efforts  have  been  devoted  to  harness  the  power  of  parallel 
computer  platforms.  While  many  CFD  and  CSD  solvers  have  been  ported  to  parallel 
machines,  grid  generators  have  lagged  behind.  For  applications  where  remeshing  is  an 
integral  part  of  simulations,  e.g.  problems  with  moving  bodies  or  changing  topologies, 
the  time  required  for  mesh  regeneration  can  easily  consume  more  than  50%  of  the  total 
time  required  to  solve  the  problem.  Faced  with  this  situation,  a  number  of  efforts  have 
been  reported  on  parallel  grid  generation. 

The  two  most  common  ways  of  generating  unstructmed  grids  are  the  Advancing  Front 
Technique  (AFT),  and  the  Generalized  Delaunay  Triangulation  (GDT).  The  AFT  in¬ 
troduces  one  element  at  a  time,  while  the  GDT  introduces  a  new  point  at  a  time.  Thus, 
both  of  these  techniques  are,  in  principle,  scalar  by  nature,  with  a  large  variation  in  the 
number  of  operations  required  to  introduce  a  new  element  or  point.  While  coding  and 
data  structures  may  influence  the  scalar  speed  of  the  ‘core’  AFT  or  GDT,  one  often 
finds  that  for  large-scale  applications,  the  evaluation  of  the  desired  element  size  and 
shape  in  space,  given  by  background  grids,  sources  or  other  means  consumes  the  largest 
fraction  of  the  total  grid  generation  time.  Unstructured  grid  generators  based  on  the 
AFT  may  be  parallelized  by  invoking  distance  arguments,  i.e.  the  introduction  of  a 
new  element  only  affects  (and  is  affected  by)  the  immediate  vicinity.  This  allows  for 
the  introduction  of  elements  in  parallel,  provided  that  sufficient  distance  lies  between 
them. 
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A  convenient  way  of  delimiting  the  possible  zones  where  elements  may  be  introduced 
by  each  processor  is  via  boxes.  These  boxes  may  be  obtained  in  a  variety  of  ways, 
i.e.  via  bins,  binary  recursive  trees,  or  octrees.  We  have  found  the  octree  to  be  the 
best  of  these  possibilities,  particularly  for  grids  with  a  large  variation  of  element  size. 
In  order  to  recover  a  parallel  gridding  procedure  that  resembles  closely  the  advancing 
front  technique  on  scalar  machines,  only  the  boxes  covering  the  active  front  in  regions 
where  the  smallest  new  elements  are  being  introduced  are  considered. 

At  the  end  of  each  parallel  gridding  pass,  each  one  of  the  boxes  gridded  can  have  an 
internal  boundary  of  faces.  For  a  large  number  of  boxes,  this  could  result  in  a  very 
large  number  of  faces  for  the  active  front.  This  problem  can  be  avoided  by  shifting  the 
boxes  slightly,  and  then  regridding  them  again  in  parallel.  This  simple  technique  has 
the  effect  of  eliminating  almost  all  of  the  faces  between  boxes  with  a  minor  modification 
of  the  basic  parallel  gridding  algorithm. 

3.  CFD-CSD-CTD  Integration 

Fluid-Structure-Thermal  Interaction  simulations  represent  one  of  many  possible  inter¬ 
disciplinary  application  areas.  One  way  to  solve  such  interdisciplinary  problems  in  a 
cost-effective  manner  is  via  the  so-called  loose  coupling  algorithm.  Each  core  discipline 
code  is  changed  as  little  as  possible,  and  the  transfer  of  information  between  regions 
or  smrfaces  is  carried  out  by  an  independent  general-purpose  library.  Results  from  a 
vast  range  of  simulations  carried  out  over  the  last  decade  indicate  that  the  proposed 
approach  offers  a  convenient  and  cost-effective  way  of  coupling  computational  fluid 
dynamics  (CFD),  computational  structural  dynamics  (CSD),  computational  thermal 
dynamics  (CTD)  codes  without  a  complete  re-write  of  them. 

The  most  important  question  is  how  to  combine  these  different  disciplines  to  arrive 
at  an  accurate,  cost-effective,  and  modular  simulation  approach  that  can  handle  an 
arbitrary  number  of  disciplines  at  the  same  time.  It  is  clear  that  any  multidisciplinary 
capability  must  have  the  ability  to  quickly  switch  between  approximation  levels,  models, 
and  ultimately  codes,  and  only  those  approaches  that  allow  a  maximum  of  flexibility, 
i.e.: 

-  Linear  and  nonlinear  CFD,  CSD  and  CTD  models; 

-  Different,  optimally  suited  discretizations  for  CFD,  CSD  and  CTD  domains; 

-  Modularity  in  CFD,  CSD  and  CTD  models  and  codes; 

-  Fsist  multidisciplinary  problem  definition;  and 

-  Fully  automatic  grid  generation  for  arbitrary  geometrical  complexity. 

We  developed  a  loose  coupling  approach,  which  includs  the  direct  coupling  in  time 
of  explicit  CFD  and  CSD  codes  and  the  incremental  load  approach  of  steady  aero- 
and  hydro-elasticity.  The  variables  on  the  boundaries  are  transferred  back  and  forth 
between  the  different  codes  by  a  master  code  that  directs  the  multi-disciplinary  run. 
Each  code  (CFD,  CSD,  CTD, ..)  is  seen  as  a  subroutine,  or  object,  that  is  called  by  the 
master  code,  or  as  a  series  of  processes  that  communicate  via  message  passing.  This 
implies  that  the  transfer  of  geometrical  and  physical  information  is  performed  between 
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the  different  codes  without  affecting  their  efficiency,  layout,  basic  functionality,  and 
coding  styles.  At  the  same  time,  different  CSD,  CTD  or  CFD  codes  may  be  replaced, 
making  this  a  very  modular  approach.  This  allows  for  a  straightforward  re-use  of  ex¬ 
isting  codes  and  the  choice  of  the  ‘best  model’  for  a  given  application.  The  information 
transfer  software  may  be  developed,  to  a  large  extent,  independently  from  the  CSD, 
CTD  and  CFD  codes  involved,  again  leading  to  modularity  and  software  reuse.  For 
this  reason,  this  approach  is  favoured  for  industrialization.  Indeed,  cosiderable  effort 
has  been  devoted  to  develop  general,  scalable  information  transfer  libraries. 

For  the  CFD-CSD-CTD  integration,  we  have  made  great  improvements  on  the  coupling 
schemes  and  information  trasfer,  including  surface  to  surface  interpolation  and  load  and 
flux  transfer,  over  the  course  of  this  research  effort. 

4.  TREATMENT  OF  MOVING  SURFACES/BODIES 

Any  fluid-structure  interaction  simulation  with  considerable  deformation  of  the  struc¬ 
ture  will  require  a  flow  solver  that  can  handle  the  arbitrary  deformation  of  surfaces 
in  time.  The  treatment  of  these  moving  surfaces  is  different  depending  on  the  mesh 
type  chosen.  For  body-conforming  grids  the  external  mesh  faces  match  up  with  the 
surface  (body  smfaces,  external  surfaces,  etc.)  of  the  domain.  This  is  not  the  case  for 
the  embedded  approach  (also  known  as  flcticious  domain,  immersed  boimdary  or 
Cartesian  method),  where  the  surface  is  placed  inside  a  large  mesh  (typically  a  regular 
parallelepiped),  with  special  treatment  of  the  elements  close  to  the  surfaces.  Consid¬ 
ering  the  general  case  of  moving  or  deforming  surfaces  with  topology  change,  both 
approaches  have  complementary  strengths  and  weaknesses: 

a)  Body-Conforming  Moving  Meshes:  the  PDEs  describing  the  flow  need  to  be  cast  in 
an  arbitrary  Lagrangean-Eulerian  (ALE)  frame  of  reference,  the  mesh  is  moved  in  such 
a  way  as  to  minimize  distortion,  if  required  the  topology  is  reconstructed,  the  mesh 
is  regenerated  and  the  solution  reinterpolated.  All  of  these  steps  have  been  optimized 
over  the  course  of  the  last  decade,  and  this  approach  has  been  used  extensively.  The 
body-conforming  solution  strategy  exhibits  the  following  shortcomings:  the  topology 
reconstruction  can  sometimes  fail  for  singular  surface  points;  there  is  no  way  to  remove 
subgrid  features  from  surfaces,  leading  to  small  elements  due  to  geometry;  reliable  par¬ 
allel  performance  beyond  16  processors  has  proven  elusive  for  most  general-purpose 
grid  generators;  the  interpolation  required  between  grids  invariably  leads  to  some  loss 
of  information;  and  there  is  an  extra  cost  associated  with  the  recalculation  of  geom¬ 
etry,  wall-distances  and  mesh  velocites  as  the  mesh  deforms.  On  the  other  hand,  the 
imposition  of  boundary  conditions  is  natural,  the  precision  of  the  solution  is  high  at 
the  boundary,  and  this  approach  still  represents  the  only  viable  solution  for  problems 
with  boundary  layers. 

b)  Embedded  Fixed  Meshes:  the  mesh  is  not  body-conforming  and  does  not  move. 
Hence,  the  PDEs  describing  the  flow  can  be  left  in  the  simpler  Eulerian  frame  of 
reference.  At  every  timestep,  the  edges  crossed  by  CSD  faces  are  identifled  and  proper 
boimdary  conditions  are  applied  in  their  vicinity.  While  used  extensively  this  solution 


9 


strategy  also  exhibits  some  shortcomings:  the  boundary,  which  has  the  most  profound 
influence  on  the  ensuing  physics,  is  also  the  place  where  the  worst  elements  are  found; 
at  the  same  time,  near  the  boundary,  the  embedding  boundary  conditions  need  to  be 
applied,  reducing  the  local  order  of  approximation  for  the  PDE;  no  stretched  elements 
can  be  introduced  to  resolve  boundary  layers;  adaptivity  is  essential  for  most  cases;  and 
there  is  an  extra  cost  associated  with  the  recalculation  of  geometry  (when  adapting) 
and  the  crossed  edge  information. 

Extensions  and  improvements  of  an  adaptive  embedded  unstructured  grid  technique 
are  made  over  past  three  years.  These  include:  improvements  in  speed  via  suitable 
data  structures,  treatment  of  multimaterial  problems,  in  particular  water/air,  the 
option  to  treat  dispersed  particles  in  the  context  of  embedded  surfaces,  a  direct  link 
to  Discrete  Particle  Methods  (DPM),  a  volume  to  surface  meshing  technique  that 
obtains  body-fitted  grids  by  post-processing  adaptive  embedded  grids;  8ind  links  to 
simplified  CSD  models.  Numerous  examples  demonstrate  the  techniques  developed. 
Prompted  by  the  need  to  solve  shock-structure  interaction  problems  and  the  inability 
of  Computational  Structural  Dynamics  (CSD)  codes  to  ensure  strict  no-penetration 
during  contact,  a  simple  embedded  technique  operating  on  adaptive,  unstructured  grids 
was  implemented.  The  essential  elements  of  this  technique  may  be  summarized  as 
follows: 

a)  The  key  modification  of  the  original,  body  fitted  edge-based  solver  was  the  removal 
of  all  geometry-parameters  (essentially  the  area  normals)  belonging  to  edges  cut  by 
embedded  surface  faces. 

b)  Several  techniques  to  improve  the  treatment  of  boundary  points  close  to  the 
immersed  surfaces  were  implemented.  Higher-order  boundary  conditions  were  also 
achieved  by  duplicating  crossed  edges  and  their  endpoints. 

c)  Geometric  resolution  and  solution  accuracy  were  enhanced  by  adaptive  mesh  re¬ 
finement  that  was  based  on  the  proximity  to  or  the  curvature  of  the  embedded  CSD 
surfaces. 

d)  In  order  to  save  work,  user-defined  or  automatic  deactivation  for  regions  inside 
immersed  solid  bodies  was  employed. 

The  technique  has  been  extremely  successful  for  complex  fluid-structme  interaction 
cases,  as  well  as  some  external  aerodynamics  cases. 

5.  OPTIMAL  SHAPE  DESIGN 

The  relentless  advance  in  numerical  methods  and  computer  power  has  made  accurate 
flow  simulations  of  realistic  geometries  a  reality.  Such  simulations  are  increasingly 
reducing  the  amount  of  lengthy  (and  costly)  experiments  in  the  aerospace,  car,  train 
and  shipbuilding  industries,  substituting  them  for  high  fidelity  CFD  runs.  This  way 
of  utilizing  CFD  is  nothing  more  than  an  exchange  of  real  for  virtual  experiment. 
However,  CFD  and  its  underlying  mathematics  offers  the  possibility  to  step  beyond 
the  capabilities  of  any  experiment.  While  the  experiment  (or  stand-alone  CFD  run) 
only  measures  the  performance  of  the  product  ‘as  is’,  numerical  methods  can  also 
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predict  the  effect  of  changes  in  the  shape  of  the  product.  This  has  led,  over  the  last 
decade,  to  a  large  body  of  literature  on  optimal  shape  design. 

A  complete  CFD  design  methodology  is  investigated.  The  main  components  of  this 
methodology  are:  a)  a  general  edge-based  compressible/incompressible  flow  solver,  b) 
a  continuous  adjoint  formulation  for  the  gradient  computations,  c)  a  steepest  descent 
technique  for  the  change  of  design  variables,  d)  evaluation  of  the  gradient  of  the  dis¬ 
cretized  flow  equations  with  respect  to  mesh  by  finite  differences;  e)  a  CAD-free  pseudo¬ 
shell  surface  parametrization,  allowing  every  point  on  the  surface  to  be  optimized  to  be 
used  as  a  design  parameter,  and  f)  a  level  type  scheme  for  the  movement  of  the  interior 
points. 

An  optimal  shape  design  procedure  based  on  the  solution  of  the  continuous  adjoint  has 
been  implemented  in  FEFLO.  Several  results,  spanning  compressible  and  incompress¬ 
ible,  viscous  and  inviscid  flow,  demonstrate  the  usefulness  of  the  capabilities  developed. 

6.  EFFICIENT  USE  OF  SUPERCOMPUTING  HARDWARE 

Despite  the  striking  successes  reported  to  date,  only  the  simplest  of  all  solvers:  explicit 
timestepping  or  implicit  iterative  schemes,  perhaps  with  multigrid  added  on,  have  been 
ported  without  major  changes  and/or  problems  to  massively  parallel  machines  with 
distributed  memory.  Many  code  options  that  are  essential  for  realistic  simulations  are 
not  easy  to  parallelize  on  this  type  of  machine.  Among  these,  we  mention  local  and 
global  remeshing,  repeated  h-refinement,  such  as  required  for  transient  problems,  con¬ 
tact  detection  and  force  evaluation,  some  preconditioners,  applications  where  particles, 
flow,  and  chemistry  interact,  fluid-structure  interaction  with  topology  change  and,  in 
general,  applications  with  rapidly  varying  load  imbalances.  Even  if  99%  of  all  opera¬ 
tions  required  by  these  codes  can  be  parallelized,  the  maximum  achievable  gain  will  be 
restricted  to  1:100.  If  we  accept  as  a  fact  that  for  most  large-scale  codes  we  may  not 
be  able  to  parallelize  more  than  99%  of  all  operations,  the  shared  memory  paradigm, 
discarded  for  a  while  as  non-scalable,  will  make  a  comeback.  It  is  far  easier  to  paral¬ 
lelize  some  of  the  more  complex  algorithms,  as  well  as  cases  with  large  load  imbalance, 
on  a  shared  memory  machine.  And  it  is  within  present  technological  reach  to  achieve 
a  100  processor,  shared  meniory  machine  (128  has  been  a  reality  since  2000). 

7.  EXAMPLES 

The  loose  coupling  methodology  has  been  applied  over  the  last  five  years  to  a  number  of 
problems.  We  include  here  some  recent  examples,  going  from  simple  rigid-body  CSD 
motion  to  highly  nonlinear,  fragmenting  (i.e.  topology-changing)  solids.  For  more 
examples,  including  validation  and  comparison  to  experiments. 

7.1  Series-60  Hull:  The  first  example  considers  the  steady  (incompressible)  flow  past  a 
typical  shiphull.  The  hull  is  allowed  to  sink  and  trim  due  to  the  fluid  forces.  The  final 
position  and  inclination  (trim)  of  the  hull  are  obtained  iteratively.  In  each  iteration,  the 
steady  flow  is  computed,  the  forces  and  moments  evaluated,  and  the  ship  repositioned. 
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The  mesh  is  typically  moved.  Should  the  need  arise,  a  local  or  global  remeshing  is 
invoked  to  removed  elements  with  negative  volumes.  The  geometry  considered  is  shown 
in  Figure  la.  The  mesh  consisted  of  approximately  400,000  elements.  Figures  lb,c 
depict  the  convergence  of  the  computed  sinkage  and  trimm  with  respect  to  the  number 
of  iterations.  Figures  ld,e  present  a  comparison  of  the  computed  sinkage  and  trimm 
with  experimental  data.  Figures  lf,g  show  a  comparison  of  the  computed  wave  drag 
coefficient  with  experimental  data  for  the  model  fixed  and  free  to  sink  and  trimm 
respectively.  A  run  of  this  kind  can  be  obtained  in  less  than  an  horn:  on  a  leading-edge 
PC. 
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Figure  lb,c  Series  60  Hull;  Convergence  of  Sinkage  and  Trim 


Figme  ld,e  Series  60  Hull:  Sinkage  and  Trim  vs.  Froude-Nr. 
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Figure  lf,g  Series  60  Hull;  Wavedrag  for  Fixed  and  Free  Model 

7.2  Nose-Cone:  Figiure  2  shows  results  for  a  proposed  nose-cone  experiment.  The 
CFD  part  of  the  problem  was  computed  using  FEFL098,  and  the  CSD  and  CTD  with 
COSMIC-NASTRAN.  Details  of  the  flow  solver  may  be  inferred  from.  The  incoming 
flow  was  set  to  Moo  =  3.0  at  an  angle  of  attack  of  a  =  10°.  The  Reynolds  number 
was  approximately  Ee  =  2  •  10®,  based  on  the  length  of  the  cone.  The  solution  was 
initiated  by  converging  the  fluid-thermal  problem,  without  any  structural  deformation. 
Thereafter,  the  fluid-structure-thermal  problem  was  solved.  Converegence  was  achieved 
after  10  cycles.  It  is  interesting  to  note  that  the  converegence  is  markedly  slower 
than  that  achieved  for  fluid-structure  (i.e.  aeroelastic)  problems.  This  is  due  to  the 
interplay  of  temperature  advection  in  the  flow  domain  and  conduction  in  the  solid, 
whose  counteracting  effects  have  to  be  balanced  out. 


Figure  2a  Nose-Cone:  Surface  Grids  for  CFD  and  CSD/CTD 


7.3  Fragmenting  Weapon:  The  third  case  considered  is  that  of  a  fragmenting  weapon. 
The  detonation  and  shock  propagation  was  modeled  using  a  JWL  equation  of  state 
with  FEFL098.  The  structural  response,  which  included  tearing  and  failme  of  ele¬ 
ments,  was  computed  using  GA-DYNA,  General  Atomics’  version  of  DYNA3D.  At  the 
beginning,  the  walls  of  the  weapon  separate  two  flow  domains:  the  inner  one,  consisting 
of  high  explosive,  and  the  outer  one,  consisting  of  air.  As  the  structure  of  the  weapon 
begins  to  fail,  fragments  are  shrunk  and  the  ensuing  gaps  are  automatically  remeshed. 
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leading  to  one  contuinuous  domain.  The  topology  reconstruction  from  the  discrete 
data  passed  to  FEFLO  from  GA-DYNA  is  completely  automatic,  requiring  no  user 
intervention  at  any  stage  of  the  simulation.  The  mesh  in  the  fluid  domain  was  adapted 
using  sources  for  geometric  fldelity  and  the  modified  H2-seminorm  error  indicator  pro¬ 
posed  in  .  The  sources  required  for  geometric  fidehty  are  constructed  automatically 
dming  the  topology  reconstruction  from  the  CSD  surface  faces.  At  the  end  of  the  run, 
the  flow  domain  contains  approximately  750  independently  flying  bodies  and  16  mil¬ 
lion  elements.  Figvue  3  shows  the  development  of  the  detonation.  The  fragmentation 
of  the  weapon  is  clearly  visible.  Figure  3c  shows  the  correlation  with  the  observed 
experimental  evidence. 


Figure  2b  Nose-Cone;  CFD/CSD/CTD  Results  Obtained 

7.4  Blast  Interaction  With  a  Generic  Ship  Hull:  Figure  4  shows  the  interaction  of  an 
explosion  with  a  generic  ship  hull.  For  this  fully  coupled  CFD/CSD  run,  the  structure 
was  modeled  with  quadrilateral  shell  elements,  the  fluid  as  a  mixture  of  high  explosive 
and  air,  and  mesh  embedding  was  employed.  The  structural  elements  were  assumed  to 
fail  once  the  average  strain  in  an  element  exceeded  60%.  As  the  shell  elements  failed, 
the  fluid  domain  underwent  topological  changes.  Figures  4a-d  show  the  structure  as 
well  as  the  pressure  contours  in  a  cut  plane  at  two  times  during  the  run.  The  influence 


of  bulkheads  on  stirface  velocity  can  clearly  be  discerned.  Note  also  the  failure  of  the 
structure,  and  the  invasion  of  high  pressure  into  the  chamber.  The  distortion  and 
inter-penetration  of  the  structiual  elements  is  such  that  the  traditional  moving  mesh 
approach  (with  topology  reconstruction,  remeshing,  ALE  formulation,  remeshing,  etc.) 
will  invariably  fail  for  this  class  of  problems. 


Figure  3a  Fragmenting  Weapon  at  310msec 


Figure  3b  Fragmenting  Weapon  at  500msec 


15 


^  70000  I 

I 


dwt  a-^Tji'a^eti  \.'clc>cltu  • 


»l 


t  i 


I  41 


1 1 1 1 


1000  1500  2000  2500  3000  350C 

llelght<gM> 


Figiire  3c  Radial  Velocity  as  a  Function  of  Fragment  Weight 


Figures  4c, d:  Surface  and  Pressure  in  Cut  Plane  at  50msec 
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8.  CONCLUSIONS  AND  OUTLOOK 

The  methodologies  and  software  required  for  Fluid-Structme- (Thermal)  interaction 
simulations  have  progressed  rapidly  over  the  last  decade.  Several  packages  (e.g. 
FEFLO/GA-DYNA3D/COSMIC-NASTRAN)  are  offering  the  possibihty  of  fully  non¬ 
linear  coupled  CSD,  CFD  and  CTD  in  a  production  environment.  Looking  towards  the 
future,  we  envision  a  multidisciplinary,  databzise-linked  framework  that  is  accessible 
from  anywhere  on  demand,  simulations  with  unprecedented  detail  and  realism  carried 
out  in  fast  succession,  virtual  meeting  spaces  where  geographically  displaced  designers 
and  engineers  discuss  and  analyze  collaboratively  new  ideas,  and  first-principles  driven 
virtual  reality. 
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